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ABSTRACT 

We describe a hybrid algorithm to calculate the formation of planets from an initial 
ensemble of planetesimals. The algorithm uses a coagulation code to treat the growth 
of planetesimals into oligarchs and explicit A'^-body calculations to follow the evolution 
of oligarchs into planets. To validate the A'^-body portion of the algorithm, we use 
a battery of tests in planetary dynamics. Several complete calculations of terrestrial 
planet formation with the hybrid code yield good agreement with previously published 
calculations. These results demonstrate that the hybrid code provides an accurate 
treatment of the evolution of planetesimals into planets. 

Subject headings: planetary systems - solar system: formation - stars: formation - 
circumstellar matter 

1. INTRODUCTION 



Rocky planet formation begins as gas and dust around a young star settle into a thin disk. 
The emergence of planets within this disk is the result of three phases of evolution (Safronov 1969; 
Weidenschilling 1980; Hayashi 1981; Wetherill & Stewart 1993; Ida & Makino 1993; Kokubo & 
Ida 1996). Initially, coagulation of dust particles causes a stochastic but steady growth in particle 
mass. The few largest bodies, or planetesimals, accumulate mass most quickly, and experience 
a runaway growth phase. As the largest objects clear out the smaller ones, the planetesimal 
growth is oligarchic, where the largest objects - "oligarchs" - become isolated from their neighbors 
and grow roughly at the same rate. In a final phase of chaotic growth, mergers of oligarchs 
lead to the formation of a few terrestrial planets around Sun-like stars well within 100 Myr (e.g., 
Weidenschilling et al. 1997; Chambers 2001; Kokubo & Ida 2002; Kokubo et al. 2006; Kominami 
& Ida 2002). 
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This scenario is bolstered by recent observations. Most, if not all, young late-type stars have a 
disk of gas and dust (Backman & Paresce 1993; Beckwith 1999; Lada 1999). At least some of these 
disks will form large planets, similar to those detected in radial velocity studies (e.g., Marcy &: 
Butler 2000). The persistence of dusty disks around older stars suggests that smaller planets form 
as well. Dust should be ejected by radiation pressure, but on-going formation of rocky planets can 
drive a collisional cascade which continually produces dust from a reservoir of small planetesimals 
(Kenyon & Bromley 2001, 2002a, 2002b). Stars with such debris disks include /? Pic (Kalas et 
al. 2000; Wahhaj et al. 2003), e Eri (Greaves et al. 1998), HR 4796A (Jayawardhana et al. 1998; 
Koerner et al. 1998) and Vega (e.g., Koerner, Sargent & Ostroff 2001). 

Interpretations of observed dusty debris disks, and our understanding of planet formation in 
general, rely heavily on numerical calculations. Two types of tools, statistical solvers and A^-body 
codes, provide complementary information about planetary disk evolution. The effectiveness of 
each tool depends on N, the number of particles they can track. Statistical methods like Safronov's 
(1969) particle-in-cell formalism work well when the particles are numerous and have small mass 
(Kokubo & Ida 1996). Current A"-body codes are unable to follow planetary growth in this 
regime. When the mass of individual objects gets large, binary interactions become important. 
Orbit evaluations then require direct A-body calculations, not the ensemble averages of particle- 
in-cell approaches. Our experience with coagulation codes (Kenyon & Luu 1998, 1999; Kenyon 
&; Bromley 2001, 2002b, 2003) suggests that we can accurately evolve objects with masses below 
10^^-10^^ g using the statistical approach. However, an A^-body code must track particles of heavier 
mass (see also Kokubo & Ida 1996, 2002). 

Massive planetesimals requiring direct A"-body evolution are relatively rare in typical planet 
formation models. This situation is fortunate because A"-body codes can not accurately track large 
numbers of particles for long periods of time. Simulations performed by Chambers (2001) and 
some runs reported here involve A ~ O(IOO) particles integrated over 100 Myr. Ida, Kokubo &: 
Kominami (2003) ran larger simulations with A = 10, 000, but only for a ~0.5 Myr time period 
and only with specialized hardware for gravitational force calculations. 

The complementary limitations of A^-body algorithms and coagulation codes call for an inte- 
gration of both methods into a single "hybrid code" (Jewell & Alexander 1996; Weidenschilling et 
al. 1997) Here we describe an algorithm that includes both a statistical component and a direct AT- 
body part. Our goal is to run self-consistent simulations of planet formation, tracking objects from 
micron-sized dust grains to Jupiter mass planets. The coagulation code, described briefly below in 
§2, accounts for gas, dust, and a swarm of lower mass planetesimals. Our A-body algorithm (§3), 
is invoked only when needed to evolve the largest planetesimals. In §4 we discuss the hybrid code 
itself, and focus on how the AT-body and coagulation parts interact. We provide tests of this code 
in §5 and give preliminary results related to the problem of terrestrial planet formation. 
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2. THE COAGULATION CODE 

The coagulation part of our hybrid code tracks dust and planetesimals in multiple annuli around 
the central star. Kenyon &; Luu (1998, 1999) and Kenyon & Bromley (2001, 2002b, 2004a) describe 
the algorithms and include a complete set of references. Here we briefly review this statistical 
algorithm. 

The physical processes which we simulate include coagulation, fragmentation, gas drag, Poynting- 
Robertson drag, and radiation pressure. We discretize the continuous distribution of particle masses 
into mass batches, assigning an integral number of particles to each batch. To allow better res- 
olution of the mass spectrum for particles evolving rapidly, adjacent batches differ in mass by a 
mass-dependent factor. We dynamically adjust the number of mass batches, as needed. Our spa- 
tial domain is cylindrical and is divided into a set of concentric annuli about the central star. The 
annuli may be set so that either they have equal width, or their boundaries are equally spaced in 
log-radius. 

Masses and particle numbers evolve according to the coagulation equations which include the 
effects of collisions, Poynting-Robertson drag, and interaction with the gas disk. Our collisions rates 
come from geometric cross sections of particles, augmented by a gravitational focusing factor for 
larger planetesimals. Collision outcomes depend on the planetesimal tensile strength, gravitational 
binding energy, and relative velocities. Depending on these parameters, collisions result in mergers, 
fragmentation, and dust production. 

We follow particle velocities statistically for each mass batch in each annulus, tracking vertical 
and horizontal velocity dispersions relative to Keplerian orbits in the central plane of the disk. 
These dispersions evolve under the influence of gas drag, Poynting-Robertson drag, viscous stirring 
and dynamical friction, according to a set of Fokker-Planck equations (Hornung, Pellat &; Barge 
1985; Wetherill & Stewart 1993; Stewart & Ida 2000; Ohtsuki, Stewart, & Ida 2002). 

To solve the coagulation and Fokker-Planck equations, we use a fourth-order Runge-Kutta 
method. The algorithm strictly conserves mass. When objects shift from one annulus to another, 
velocity dispersions are updated to conserve kinetic energy explicitly. At the inner and outer 
boundaries of the spatial domain, the mass batches reflect a steady, radial flow in the debris 
disk. The code does not enforce angular momentum conservation, nonetheless, in test runs of 10^ 
timesteps, angular momentum is conserved to better than 1%. 

The coagulation code evolves the swarm of 0(10^^) planetesimals accurately until the largest 
objects reach roughly the mass of Pluto. Because these massive objects are rare, the statistical 
model provides a poor estimate of their behavior. Although the coagulation code can still follow 
the evolution of the swarm, the orbits of the largest objects require explicit calculations using an 
iV-body code, which we describe in the next section. 
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3. THE A^-BODY CODE 

Planetary A^-body codes must be fast and accurate. Suitable solvers for the equations of 
motion are numerous, with published descriptions of symplectic integrators (Wisdom &; Holman 
1991; Kinoshita, Yoshida k. Nakai 1991; Saha & Tremaine 1992), and more general time-symmetric 
integrators (e.g., Quinlan & Tremaine 1990, Hut, Makino & McMillan 1995). Some integrators 
(Vasilev 1982; Wisdom & Holman 1991; Ida & Makino 1992; Fukushima 1996; Shefer 2002) speed 
up the calculations by tracking only deviations from a purely Keplerian orbit about the central 
star, following the proposal of Encke (1852). Other improvements include adaptive timestepping, 
allowing integrators to take large timesteps when forces are small, and to expend more computa- 
tional resources only when particles experience strong forces (Aarseth 1985; McMillan & Aarseth 
1993; Skeel & Biesiadecki 1994; Kokubo, Yoshinaga & Makino 1997; Duncan, Levison & Lee 1998). 

Our algorithm is based on the Encke method. We work with non-inertial, Keplerian frames 
about the central star and integrate the equations of motion in rectilinear coordinates defined 
relative to the Keplerian frames. These coordinates have a fixed orientation, aligned with the 
incrtial frame of the central star. We use an adaptive block timestepping scheme (e.g., McMillan & 
Aarseth 1993) and a sixth- or eighth-order accurate integrator based on Richardson extrapolation. 
We give details of the integrator below. For now we briefly outline the timestepping scheme. 

At the beginning of a timestep, each particle's accelerated reference frame is initially set to 
its instantaneous Keplerian orbit. If particles are in tight groups, their reference frames can be set 
to the Keplerian orbit of their mutual center of mass. A friends-of-friends algorithm (Huchra & 
Geller 1982, Geller k, Huchra 1983) finds any such groups with a linking parameter based on the 
Hill radii of the particles. We then integrate the equations of motion in terms of spatial variables 
in the accelerated reference frames over a single timestep of length At. Next we calculate higher 
resolution orbits by dividing the timestep into m equal substeps. In practice, m = 3 seems to 
work best. A comparison between the lower and higher resolution orbits produces a check of 
integrator convergence. We have implemented several criteria for convergence, but have found that 
an individual particle's total energy, taken as a fraction e, of its Keplerian orbital energy, gives an 
effective comparison under most circumstances. If a particle's orbit has not converged, we simply 
repeat the process at even higher time resolution. 

The strategy of a block timestepping scheme is to track individual orbits and interpolate the 
converged orbits if needed to evaluate forces on particles whose orbits have not converged. This 
procedure reduces the computational load of force calculations if the number of particles is large 
{N S> 10). For a smaller number of particles we prefer to integrate all orbits until every orbit 
has converged. Energy conservation is better in this case and there is no significant penalty in 
computational load. 

The time integrator is an ordinary differential equation (ODE) solver based on Richardson 
extrapolation. We start with a low-resolution estimate of the particle orbits over a time interval 
At using a leapfrog intergrator. Specifically, the position x and velocity w of a particle at the end 
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of a time-forward step in the leapfrog case is 



v{t + At/2) = v{t) + At/2a{t) 

x{t + At) = x{t) + Atv{t + At/2) 

v{t + At) = v{t + At/2) + At/2a{t + At) 



(1) 



where a is the acceleration, which in our case depends only on position. We label the position at 
the end of this single leapfrog timestep as xq, and note that it contains errors of the order of At^. 
Next wc divide the time interval into two equal parts and take two successive leapfrog steps to 
derive a better resolved orbital position, xi. We can further subdivide the time interval to obtain 
a position Xj, derived from 2* successive leapfrog timesteps. Our final position at the end of the 
time interval At is a linear combination of these results. 



Fourth-, sixth-, and eighth-order methods have m= 1, 2 and 3, respectively and coefficients 



We calculate the gravitational forces between particles directly. At the lowest resolution 
timesteps, we perform the 0{N'^) force evaluations. However, if we use the block timestepping 
scheme then only a few particles depend on 0{N) interparticle forces at high temporal resolution. 
To simulate a larger number of bodies, we have an 0{N\ogN) Barnes & Hut (1986) treecode, as 
described in Barton, Bromley & Geller (1999). Preliminary tests suggest that this code becomes 
competitive with a direct method only when N is greater than O(IO^). An alternative would be to 
use specialized hardware (e.g.. Hut &: Makino 1999). Here, our particle numbers are small and we 
work exclusively with the direct force solver. 

For realistic planet formation models, an A^-body algorithm must identify merger events. Wc 
check for mergers in a fast way, by assuming that any merger event which occurs during a single 
timestep can involve only two particles. For each particle we save an index number corresponding 
to some other particle which came the closest during the force calculations. This indexing reduces 
the number of calculations from 0{N'^) to 0{N). For each pair of bodies, with indices i and j, we 
check the relative radial velocity using the quantity Xij ■ Vij, where the vectors are relative position 
and velocity respectively. If the relative radial velocity is positive at the beginning of a timestep, 
then we assume no merger takes place. Otherwise, we find the minimum separation of the pair 
during the timestep, interpolating if necessary. Once we have the minimum pair separation Sij and 
relative speed Vij, a merger event is identified if 



m 




(2) 



(co,Ci) = (4/3,-1/3) 
(co,ci,C2) = (1/45,-4/9,64/45) 
(co,ci,C2,C3) = (-1/2835,4/135,-64/135,4096/2835) 



4^1 —order 
6*'*— order 
8*'*— order. 



(3) 



Sij < min Ri + Rj, — {vE,iRi + vejRj) 



(4) 
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where VE,i is the escape velocity at distance Sij from the i particle, and Ri is the particle's 
physical radius. Thus the merger cross-section is at least the physical cross-section, and is larger if 
the relative velocity is small compared to the escape velocity of either particle. 



With our coagulation code and A^-body algorithm, we track both the numerous low-mass 
planetesimals and the relatively rare high-mass bodies in a planetary disk. The interactions between 
these two populations generally cause the larger bodies to circularize in their orbits, while the 
smaller bodies tend to get gravitationally stirred, as reflected in an overall increase in eccentricity 
and inclination. Here we describe how the coagulation and iV-body components work together to 
simulate these effects. 

The multiannulus coagulation code evolves all low-mass bodies. The model grid contains 
N concentric annuli with widths 5a j centered at heliocentric distances Uj. Each annulus contains 
n{mji,t) objects of mass niji with orbital eccentricity eji{t) and inclination iji{t) in M mass batches. 
When an object in the coagulation code reaches a preset mass, it is 'promoted' into the A^-body 
code. To set the initial orbit for this object, we use the three coagulation coordinates, a, e, 
and i, and select random values for the longitude of periastron and the argument of perihelion. 
Because the annuli have finite width Saj, we set the semimajor axis of the promoted object to Up = 
aj + (0.5 — x)5aj, where x is a random number between and 1. When two or more objects within 
an annulus arc promoted to the A^-body code during the same timestcp, we restrict the choices of 
the orbital elements to minimize orbital interactions between the newly promoted A'^-bodies. 

The coagulation code also determines the simulation timestep. This interval, AT, is generally 
larger than the lowest resolution timestep of the A^-body code, so we take multiple low-resolution 
substeps, At, within that interval. The length of a substcp is set by the orbital speed of the fastest 
moving body; for nearly circular orbits this limit is roughly 1/60*^ of the orbital period. 

To calculate the effects of the swarm of low-mass bodies on the A'^-bodies, we use particle-in- 
a-box estimates. For an N-body with index j, mass rrij, eccentricity ej, inclination ij, horizontal 
velocity hj, and vertical velocity Vj, we derive the time-evolution of the orbital eccentricity and 
inclination from the Fokker-Planck formulae for the derivatives of the horizontal and vertical ve- 
locities: 



4. 



THE HYBRID CODE 





k=l m=l 



N 



M 
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for viscous stirring in the high velocity hmit and 

-^f^ = E E l-^/.'i^ {rUkmhlra - mjh]) HeiPkm) (7) 
fe=l m=l 

^^2 k=N M 

k=l m=l 

for dynamical friction in the high velocity limit. In these expressions, = {ij + i'l„J/{e'j + e|.^„), 
C = 0.2767 C^prn/ih^j + hl^f/"^, and the functions He, H-„ Jg, and are definite integrals 
(Stewart &: Ida 2000; Ohtsuki, Stewart, & Ida 2002). The subscript k is the annular index while 
subscript m is the index of a mass batch. The overlap fraction /j^ is the fraction of bodies in 
annulus A; that approach within 2.4 Rh of the iV-body. We set pm = Mj^/Vm, where is the 
total mass of bodies with rUm in annulus I and Vi is the volume of annulus I. Following Stewart & 
Ida (2000), we also set Aa = ln(A2 + 1), where 

A = 0.1886 {h + 1.25?;) v/v^^^^. , (9) 

with h? = h'] + hl^, = v] + vl^, VHjk = rnjkVKjk, r]j -^. = {rrij + mkm)/MQ, and VKjk = 
0-5{vK,j + VK,k)', the subscript K denotes a Keplerian velocity (see also Kenyon &; Bromley 2002b; 
Ohtsuki, Stewart, & Ida 2002). 

When the relative velocities of particles approach the Hill velocity, VHjk, we set 

ln(10A^(e^)V^ + l) 
C2 - ioA2(e2)i/2 

ln(10A^ + 1) 

= 10A2 (^2) 

and use 

ju2 k=Nl=M 

= E E 73/,.Ci4,, (13) 
fc=i i=i 



2 k=Nl=M 
vs,j,low 



dt 

k=l 1=1 



Y: E /.^C72(4(z2) + 0.2{er/'{^'y/'K^^, (14) 
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for viscous stirring in the low velocity limit and 



77 2 k=N l=M 



dt 



k=l 1=1 

k=Nl=M 



dt 



(15) 



(16) 



k=i 1=1 



for dynamical friction in the low velocity limit (Ohtsuki, Stewart, &; Ida 2002). In these expressions, 
C4 = 0M25pmHki{aj + ak)V^j^./{Mj + Mkm)^, Hjk is the mutual scale height, (e^) = {ej + 
eljjj/rjf jf., and (z^) = {ij + iljn)/fH,jk- '^^^ combined velocity stirring is then 



dt 



"'"'vs,j,hi9h _|_ ^"'vs,j,low 



dt 



dt 



(17) 



dt 



dt 



+ 



dt 



(18) 



for viscous stirring and 



^u2 ju2 ]u2 



dt 



dt 



+ 



dt 



(19) 



for dynamical friction. 



dt 



dt 



dt 



(20) 



To calculate the accretion rate of planetesimals onto the AT-bodies, we use the standard coag- 
ulation equation: 



drrij 
~df 



k=N l=M 



(21) 



k=l 1=1 



where Ajj^k is the normalized cross-section (see Kenyon Sz Luu 1999). 

To calculate the stirring of planetesimals by iV-bodies, we calculate the appropriate Fokker- 
Planck terms for viscous stirring and dynamical friction and add these results to the long distance 
stirring of Weidenschilling (1989) 
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dhl 



j=N, 



n 



Aa 



dt 



(22) 



where is the number of A^^-bodies and Aa = aj — Ofc. 

At the end of each coagulation timestep, we pass the stirring and accretion rates for each 
iV-body - dej/dt, dij/dt, and druj/dt - to the AT-body code. At the end of every low resolution 
timestep in the AT-body code, we modify each particle's orbit and mass to reflect these changes. 
When possible, we simply redirect the velocity vector so that the eccentricity varies independently 
of inclination and semimajor axis. The iV-body functions return an updated particle list to the 
coagulation code. New orbital positions and masses of the large objects are then reinserted into 
the coagulation grid. The coagulation calculations proceed in the standard fashion, except that the 
orbital velocities of the large-mass objects in the grid are not evolved. With a complete circuit from 
the coagulation code to the AT-body code to the coagulation code, A?^-bodies influence the evolution 
of the swarm and the swarm influences the evolution of the A/"-bodies. 



We have published elsewhere results on the performance of the coagulation code, as described 
in §2. However, the A^-body and hybrid codes are new, and our purpose here is to validate them. 
In this section we first test the A^-body code for stability, dynamic range, accuracy, and merger 
resolution, mostly following Duncan, Levison & Lee (1998) in the validation of their SyMBA algo- 
rithm. We then test the hybrid code against several simulations of terrestrial planet formation at 
1 AU (Weidenschilling et al. 1997; Chambers 2001). 

The following subsections are organized according to the type of calculation. We test stability 
during long term orbit integrations of the major planets and a "scaled outer solar system" where 
the masses of the major planets are increased by a factor of 50 (Duncan, Levison &; Lee 1998). 
Then we test dynamic range with two binary planet configurations. The code's accuracy is also 
established with a test to resolve a critical orbital separation between two planetesimals which 
determines whether their orbits will cross. In considering mergers, we reproduce the Greenzweig 
&; Lissauer (1990) results for planetary accretion rates. We also derive the integration accuracy 
required to track the collision between a massive object and a small projectile on opposing circular 
orbits at 1 AU. Finally, we simulate terrestrial planet formation, following Weidenschilling et al. 
(1997), who consider the evolution of km-sized planetesimals, and Chambers (2001), who models 
the coUisional evolution of lunar mass bodies. 



5. 



TESTS 
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Long-term Orbit Integration 

To evaluate the stability of our adaptive integrator, we evolve the four major planets in orbit 
about a stationary Sun. Figure 1 shows the behavior of the outer planets' eccentricities over a 
10 Myr period. We use two different integrators, a sixth-order accurate symplectic integrator 
(Yoshida 1990) and our adaptive code, also with a sixth-order accurate ODE solver for comparison. 
We run both codes at low and high resolution. The run times of the two codes are comparable 
at each resolution. As illustrated in the figure, the low resolution adaptive code obtains a better 
measure of eccentricity than its symplectic counterpart. One reason is that the adaptive code is 
not subject to accumulation of phase errors as in a symplectic algorithm. On the other hand, the 
adaptive code produces a secular drift in energy over the course of the simulation, whereas the 
symplectic code does not. The fractional error in total energy at 10 Myr is 10~^ and 10~^ for the 
low and high resolution adaptive code, respectively. 

Scaled Outer Solar System 

If the masses of each of the four major planets increase by a factor of 50, the outer solar system 
is unstable (sec, for example, Duncan &; Lissauer 1998). To verify that our code yields the correct 
timescale for this instability, we made tests similar to one performed by Duncan, Levison & Lee 
(1998). As expected, the semimajor axis of Jupiter's orbit shrinks by a modest fraction of the 
original semimajor axis, while Saturn is ejected on a timescale of 1,000 years. As Saturn leaves the 
Solar System, its orbit crosses and excites the orbits of Uranus and Neptune. An increase in the 
orbital eccentricity of Uranus or Neptune leads to a second ejection, as seen in Figure 2. 

Binaries 

Duncan, Levison &: Lee (1998) note that bound binary planets whose center of mass revolves 
around a star provide a strong test of the ability of an A^-body integrator to handle a wide range of 
dynamical timescales. They set up a pair of Jupiter-mass planets with a center of mass on a circular 
orbit at 1 AU. The planets orbit each other with an initial binary semimajor axis of 0.0125 AU and 
binary eccentricity of 0.6. Their simulation time was 100 yr, covering about 3,200 binary orbits. 

We reproduce the Duncan, Levison & Lee (1998) test. We can limit total energy errors to 
within a part in 10 million, while keeping the wall time well below a minute on a 1.7 GHz AMD 
Athlon processor. Figure 3 illustrates changes in the binary semimajor axis. 

The figure also demonstrates a test of our block timestepper. If a small particle is in a tightly 
bound orbit around a much more massive object, then the position of the massive particle is 
interpolated in the force calculations for the small body. Our test simulates the orbit of a massless 
"spy satellite" on a polar orbit about the Earth. Figure 3 shows that the distance of the satellite 
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from the center of the Earth varies by about 100 meters during an orbit. There is no significant 
drift in its mean altitude over the course of 100 years, corresponding to roughly 600,000 satellite 
orbits of the Earth. 

Planetesimals 

We next test the ability of the A^-body code to handle the evolution of planetesimals. First, 
we consider an isolated pair of planetesimals, with masses of 2 x 10^^ g. These bodies start on 
corotating, circular orbits near 1 AU. We specify their orbital separation in terms of their mutual 
Hill radius, Rh = [{mi + m2)/3Mo]^/^. There is a critical separation, Aacrit = 2\/3Rh, which 
determines the onset of chaotic behavior: If the planetesimals' initial orbital separation is smaller 
than Aacrit, their orbits will cross, otherwise they will simply scatter to larger orbital separations. 
Figure 4 shows the evolution of semimajor axes in cases where orbital separations are initially 
smaller than Aocrit (top panel) and greater than Aocrit (middle panel). As expected, the more 
closely-spaced pair experiences several orbit crossings while the more distant pair does not. Multiple 
calculations demonstrate that our code can resolve the critical separation Aocrit to within a few 
percent. 

Figure 4 also shows results for eight planetesimals, each with the same mass as in the preceding 
case, and initially on circular, Keplerian orbits. The separation between nearest neighbors is 
1.05Aacrit- The orbit crossings proceed from larger radius inward, as illustrated anecdotally in the 
figure (see also Weidenschilling et al. 1997; Kenyon & Bromley 2001). 

Next we use the A^'-body code to model gravitational stirring by two 2 x 10^^ g objects in a 
uniform disk of 805 lower mass, 2 x 10^4 g objects (Kokubo & Ida 1995; Weidenschilling et al. 1997; 
Kenyon &; Bromley 2001). The disk is centered at 1 AU, and is 35 Rh in annular extent, where Rh 
is the mutual Hill radius between large and small objects. To speed up the code's convergence with 
this relatively large number of particles, we soften the interaction potential on a scale comparable 
to the physical radius of the large bodies, assuming their density is 1 g cm~^. 

Figure 5 illustrates the evolution of the radial profile of (e^ + 1^)^^"^ for the case where the two 
massive bodies are separated by IORh, and are each at an orbital distance of 5Rh away from 1 AU. 
The initial eccentricity and inclinations are zero for the two massive objects, and small (i ~ 10~^, 
e ~ 10^'^) for the planetesimals. The A^-body results here compare well with those from the pure 
coagulation code (Kenyon & Bromley 2001), except that the A^-body calculations produce more 
structure in the profile, particularly an enchanced scattering of the swarm away from the more 
massive bodies. This difference, a broadening of the profile as compared to the coagulation results, 
is expected from previous AT-body calculations (Kokubo &; Ida 1995). 
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Mergers 

One test of our merger algorithm is provided by Greenzweig k. Lissauer (1990), who studied 
gravitational focusing by a single planet in a field of test particles, all orbiting a 1 M© star. We 
consider one of their configurations, a 10~^ planet on a circular orbit at 1 AU and a set of 

test particles on orbits with eccentricity e = 0.007, inclination i = 0.2 degrees, and semimajor axes 
distributed in two rings between 0.977 AU and 1.023 AU. We evolve the system through a single 
close encounter between the planet and each test particle, to determine the fraction of test particles 
accreted by the planet. 

When the planet has a physical radius of 10^ km, Greenzweig & Lissauer (1990) report an 
accretion fraction /acc = 0.140. In their validation of the SyMBA code, Duncan, Levison &; Lee 
(1998) estimate /acc = 0.10 ± 0.02. For a planetary radius of 5,200 km, the derived accretion 
fractions are /acc = 0.009 (Greenzweig & Lissauer 1990) and /acc = 0.008 it 0.002 (Duncan, Levison 
& Lee 1998). Using the merger criterion specified in §3 (eq. [4]), our results are /acc = 0.138ib0.002 
(Poisson errors with 50,000 test particles) for the 10^ km radius planet and /acc = 0.008 it 0.001 for 
the 5,200 km radius planet. 

A more extreme test of the merger algorithm is to simulate the high-speed collision between 
two counterrotating objects. We consider a small (cm-size) projectile and a larger target of some 
specified radius, r, both with a density of 1 g/cm'^. The two bodies initially are at 1 AU on opposite 
sides of the Sun and are on colliding circular orbits. The code calculates trajectories with a fixed 
number of timesteps over a time interval which is randomly distributed at values just greater than 
0.25 yr. We run repeated trials and find the minimum number of low-resolution timesteps, n^s, 
required to achieve a 95% success rate in detecting collisions. 

Figure 6 illustrates our results. For targets with r < 1, 000 km, 71^95 is determined by the errors 
in the interpolation of position between the endpoints of each timestep, and scales approximately as 
J^tQS ~ r~^/'^. In our planet formation simulations reported here, we typically take O(IOO) timesteps 
per orbit and therefore can expect to resolve high-speed collisions between 1,000 km objects. If 
gravitational interactions between target and projectile become important, as in cases with large 
bodies and slow impact speeds, then the adaptive integrator takes high-resolution timesteps, in- 
creasing the accuracy of the merger algorithm. This effect is illustrated in Figure 6, which shows 
a dramatic decrease in njgs for targets of size 100 km or larger, depending on the adaptive code's 
error tolerance. In practice, our iV-bodies typically have radii greater than 1,000 km and have im- 
pact speeds which are more than an order of magnitude smaller than in a counterrotating collision. 
Thus, our code should accurately resolve mergers in planet formation simulations. Next, we test 
this assertion explicitly. 
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Terrestrial Planet Formation 

To illustrate the behavior of the hybrid code for less idealized conditions, we consider calcu- 
lations of terrestrial planet formation at 1 AU. We compare our results with two sets of published 
calculations. Spaute et al. (1991) and Weidenschilling et al. (1997) use a multi-annulus code which 
evolves the masses and orbital properties of planetesimals with a coagulation code and follows the 
evolution of discrete planetary embryos with a Monte Carlo algorithm. As in our calculations, their 
code allows interactions between the planetesimals and the embryos. They evolve planetesimals 
for 1 Myr in two configurations, 0.86-1.14 AU and 0.5-1.5 AU. Here, we use a grid at 0.84-1.16 
AU and compare the Weidenschilling et al. (1997) results with our calculations from the pure 
coagulation code and from the hybrid code. 

Chambers (2001) uses a three-dimensional AT-body code to consider the final phase of terres- 
trial planet formation at 0.4-2 AU. These calculations begin with 150-160 lunar-mass to Mars-mass 
embryos and follow the collisional evolution for 200-500 Myr. Here, we adopt a grid at 0.4-2 AU 
and compare the outcomes of two pure A^-body calculations and many hybrid calculations with the 
Chambers (2001) results. 

We begin with results for a = 0.84-1.16 AU. In these calculations, the grid has 32 radial zones 
each containing an initial distribution of planetesimals with radii of 1-4 km. The planetesimals 
have initial surface density E = 16(a/l AU)~^/^ g cm~^, eccentricity cq = 10~^, and inclination iq 
= 6 X 10~^. The calculations include gas drag but do not allow fragmentation. 

The initial stages of pure coagulation and hybrid calculations are identical. In a few thousand 
years, objects grow to radii of 10-20 km. After ~ 10'^ yr, objects reach sizes of 100-300 km. Because 
the timescales for dynamical friction and viscous stirring are shorter than the growth timescales, 
the largest objects have nearly circular orbits while the smallest objects have eccentric orbits. 
Runaway growth then yields a handful of 1000 3000 km objects. As runaway growth proceeds, 
viscous stirring continues to heat up the orbits of the smallest objects faster than the smallest 
objects damp the velocities of the largest objects. Thus, gravitational focusing factors decrease, 
accretion slows, and runaway growth ends. The evolution then enters a period of 'oligarchic growth,' 
where all of the largest objects grow at roughly the same rate (Kokubo k, Ida 1996,1998,2002). 

During oligarchic growth, the growth of the largest objects in the hybrid calculations differs 
from the path followed in the pure coagulation models. In the hybrid models, oligarchs grow slowly 
until their orbits begin to cross (Kokubo & Ida 2002; Kominami k, Ida 2002; Kenyon & Bromley 
2006). Once orbits cross, chaotic growth leads to a rapid merger rate and the formation of several 
'super-oligarchs' that accrete most of the leftover planetesimals. The super-oligarchs accrete some 
of the remaining lower mass oligarchs and scatter the rest out of the grid (Kenyon &; Bromley 
2006). 

In pure coagulation models, the orbits of oligarchs arc not allowed to cross. These oligarchs 
slowly accrete all of the planetesimals within their gravitational reach. Oligarchs rarely accrete 
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other oligarchs. After 10-100 Myr, pure coagulation models have more lower mass oligarchs than 
the hybrid models. However, the largest oligarchs are less massive than their counterparts in the 
hybrid calculations. 

Figures 7 and 8 compare some results from two calculations at specific times. Figure 8 of 
Wcidcnschilling et al. (1997) shows results for a similar calculation. In every case, the mass of the 
largest object is ~ 5 — 6 x 10^^ g at 10^ yr, ~ 3 — 4 x 10^^ g at 10^ yr, and ~ 10^^ g at 10® yr. 
At early times, e and i decline monotonically with increasing mass. At later times, e and i have 
bimodal distributions: objects with m < 10^^ g have large e and i, while more massive objects have 
much smaller e and i. The orbital anisotropy, i/e, follows a similar evolution, with z/e ~ 0.4 at 
early times and i/e ^ 0.1-0.15 at late times. Near the end of the calculations, the largest objects 
lie in a flattened disk, with i/e ^ 0.02-0.03 (see also Weidenschilling et al. 1997). 

Figure 9 illustrates the evolution of the eccentricity in more detail. At 10^ yr, dynamical 
friction maintains low e for the largest bodies. Because planets form fastest at the inner edge of 
the grid, small objects at the inner edge have larger e than small objects at the outer outer edge. 
At 10^ yr, dynamical friction still maintains low e for the largest bodies. Because most of the small 
objects have grown to 100 km radius, these objects have slightly smaller e than other objects with 
m < 10^^ g. At 10® yr, objects with m > 10^® g have low e ~ 0.02, while lower mass objects have 
higher e ~ 0.1. These results compare favorably with Figure 9 of Weidenschilling et al. (1997). 

We begin our version of the Chambers (2001) calculations in a grid of 40 annuli at a = 0.4-2 
AU. For the hybrid calculations, the anmili contain an initial distribution of planctesimals with radii 
of 4-15 km. The pure A'^-body model starts with 160 'moons' with a mass of 10^® g. In both cases, 
the objects have cq = 10~^ and io = 6 x 10~^. To provide some contrast with previous calculations, 
we adopt an initial surface density of planctesimals or moons, S = x'Eo{a/l AU)~^ with Eq = 8 g 
cm~^ and x = 0.25-2, instead of the usual S oc a""^/^. Although the final configuration of planets 
depends on the initial surface density gradient, the general evolution is fairly independent of S. 

The large radial extent of these calculations allows us to illustrate the sensitivity of planet 
formation to the heliocentric distance. Because the timescale for growth by coagulation is t oc 
E/P oc a^/^, planets grow fastest at the inner edge of the grid (Lissaucr 1987). In hybrid models 
with X = 1, objects with radii of ^ 200 km form in ~ 10'^ yr at 0.4 AU, in ~ 10^ yr at 0.95 AU, 
and in ~ 10^ yr at 2 AU. Oligarchs with masses of ~ 10^® g form on a timescale 

to ~ 1.5 x 10^x-^{a/l AU)5/2 . (23) 

Once large objects start to form, the transition from runaway to oligarchic to chaotic growth 
proceeds in several waves propagating from the inner disk to the outer disk (Kenyon & Bromley 
2006). As oligarchs start to form at the outer edge of the grid, dynamical interactions between 
oligarchs begin at the inner edge. Several orbit crossings lead to dynamical interactions between 
all oligarchs and a rapid increase in the merger rate. A few large oligarchs gradually accrete many 
of the smaller oligarchs, leading to a configuration with several Earth-mass planets and a few 
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Mars-mass 'leftovers.' 

In the moons calculations, dynamical interactions between the 160 original oligarchs dominate 
the entire evolutionary sequence. In the first ~ 10^ yr, five mergers start to produce large oligarchs 
in the inner disk. The number of oligarchs declines to 146 in 10^ yr, to 115 in 1 Myr, and to 55 in 
10 Myr. By 10 Myr, the largest objects in the inner disk reach masses of 0.2-0.3 and slowly 
accrete the remaining moons. In the outer disk, the smaller surface density and viscous stirring 
lead to a smaller merger rate. After ~ 100 Myr, only ^ 20 oligarchs remain in the outer disk. 

To compare our results with Chambers (2001), we consider measures of the orbital elements 
of final 'solar systems' from several simulations. In our direct A^-body and hybrid calculations, the 
time variation of the mass-weighted eccentricity e (Fig. 10) follows a standard pattern (see Fig. 5 of 
(Chambers 2001)). During runaway and oligarchic growth, protoplanets stir their surroundings and 
e increases. At 10-100 Myr, the orbits of oligarchs cross, which further excites orbital eccentricity 
and produces sharp peaks in e. Mergers between oligarchs reduce e. Once a few remaining planets 
have fairly stable orbits, e settles to 0.05-0.15. 

Hybrid models yield larger variations in e than A-body calculations of large objects. Because 
objects grow to the promotion mass throughout a hybrid calculation, the merger phase lasts longer 
and produces more frequent resonant interactions at later times compared to pure A-body models. 
Thus, we often observe several peaks in the evolution of e for hybrid models compared to a single 
peak in pure AT-body calculations. 

Figure 11 shows the relation between the eccentricity and mass of the planets produced in our 
calculations. Our models yield a range in mass from ~ 0.01 to ~ 2-3 M®. For ~ 30 planets 
with masses, m > 0.1 M^, the orbital eccentricity is not correlated with the mass. Chambers 
(2001) derived a similar result for ~ 50 planets. However, our calculations yield a reasonably large 
ensemble of lower mass planets with m ~ 0.01-0.1 M®. In the full ensemble of planets with m ~ 
0.01-3 M®, there is a small, but significant trend of decreasing eccentricity with increasing planet 
mass. From the Spearman rank and Kendall's r tests (Press et al. 1992), the probability that the 
derived distribution of m{e) is random is ~ 2 x 10~^. 

To compare this result to our Solar System, we restrict the mass range to 0.05-2 M^. For our 
ensemble of ~ 40 planets, the Spearman rank and Kendall's r probabilities are ~ 0.1. These tests 
yield probabilities of ~ 0.15-0.2 for the 4 terrestrial planets in our Solar System. We conclude that 
the weak correlation of e with m is consistent with architecture of our Solar System. 

As a final comparison with Chambers (2001), we calculate several statistical parameters to 
characterize the final configurations of our models. For a system with N oligarchs, 

Sm = mi/mt (24) 

measures the fraction of the mass in the largest object, where mi is the mass of the largest oligarch 
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and mt is the total mass in the grid. The orbital spacing statistic is 



6 / a. 




'2M 



1/4 



Ss = 



N-1 \a, 



(25) 



where fh is the average mass of an oligarch, nic is the mass of the central star, a^in is the minimum 
semimajor axis of an oligarch, and Umax is the maximum semimajor axis of an oligarch. To measure 
the degree of orbital excitation of a planet, Laskar (1997) describes the angular momentum deficit, 
the difference between the z-component of the angular momentum of an orbit and a circular orbit 
with the same semimajor axis. Chambers (2001) generalizes this definition as the sum over all 
planets. 



measures whether mass is concentrated in a few massive objects (large Sc] as in the Earth and 
Venus) or many low mass objects (small Sc, as in the Kuiper belt and scattered disk of the outer 
solar system). 

Table 1 compares results for our 14 hybrid calculations and our two 'moons' calculations with 
statistics for the inner Solar System. For calculations with x = 1, our moons calculations yield 
roughly the same outcome as Chambers (2001). The mass of the largest object and the statistics 
agree well with those in Table 1 of Chambers (2001). Because we began with a shallower surface 
density gradient, our planetary systems have more massive objects and more oligarchs at larger 
heliocentric distances than Chambers (2001). 

Our hybrid models for x = 1 yield results similar to the moons calculations. Our most massive 
object has a mass of 1.1-1.8 M® and contains 40% to 80% of the initial mass in the grid. The 
calculations typically produce 1-2 other massive objects comparable in mass to Venus and several 
less massive objects with masses more similar to Mars. For models with longer evolution times, 
the number of Mars-mass objects declines considerably from 100 Myr to 500 Myr. The statistical 
parameters agree well with the Chambers (2001) results. The planets in all of our model solar 
systems have more eccentric orbits than the planets in the Solar System. However, the other 
statistics agree rather well with those in the Solar System, including one model that is more 
concentrated than the Solar System. 

The outcomes of our calculations are sensitive to the initial mass in the grid. Models with x 
= 0.25 produce very low mass planets, ~ 0.1 M^, compared to models with x = 1-2, ~ 1-2 M®. 
The lower mass models also yield more planets on more circular orbits and have a smaller fraction 
of the total mass in the largest object. 



Sd = 




(26) 



Finally, the mass concentration statistic 




(27) 
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6. CONCLUSION 

In this paper we describe a hybrid A^-body coagulation code for planet formation. We put the 
iV-body part through a battery of tests to assess its performance in planetary dynamics simulations. 
Validation tests of the coagulation algorithm appear elsewhere (e.g., Kenyon & Bromley 2001). We 
provide tests of the hybrid method by comparing it to coagulation simulations and A'^-body output 
separately. We demonstrate that all three methods quantitatively reproduce published results for 
the formation of rocky planets by mergers of 0(10"^^) km-sized (Weidenschilling et al. 1997) or 
O(IOO) moon-sized (Chambers 2001) planetesimals. Table 2 gives a summary of these tests, with 
references to other published work. 

Although our intent here is simply to describe our method, the hybrid code now has the 
potential to match observations of the Solar System and debris disks (e.g. Kenyon & Bromley 
2004a,b). In addition to results on terrestrial planets in §5, Kenyon & Bromley (2006) summarizes 
results on the transition from oligarchic to chaotic growth in the terrestrial zone. Kenyon & 
Bromley (2004c) shows how a stellar encounter can produce Scdna-like orbits in the outer Solar 
System and outlines how observations might distinguish between locally-produced and captured 
Sednas. Future papers will describe more complete results on planet formation in the terrestrial 
zone and the trans-Neptunian region. 

We acknowledge a generous allotment, ~ 20 cpu years, of computer time on the Silicon Graphics 
Origin-2000 'Alhena' through funding from the NASA Offices of Mission to Planet Earth, Aero- 
nautics, and Space Science. Advice and comments from M. Geller are greatly appreciated! The 
NASA Astrophysics Theory Program supported part of this project through grant NAG5-13278. 
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Table 1. Statistics for Planetary Systems at 0.4-2 AU 





U„oi (Myr) 
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Sd 


Sc 


0.25 
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22.4 


1.00 
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1.00 


500 


2 


2.1 


0.778 


36.1 


0.0228 


119.4 


2.00 
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2.00 
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Table 2. Summary of Numerical Tests 



Test name 


Result (this paper) 


Result (other references) 


Tesfti of in.fcyrdtion (iccuracy: 
Giant planet orbit integration 
Scaled Outer Solar System 


Figure 1 
Figure 2 


Figure 5, Duncan, Levison & Lee (1998) 


Tests of dynamic range: 
Binary Jupiters 
Earth satellite 


Figure 3 
Figure 3 


Figure 4, Duncan, Levison & Lee (1998) 


Planetesimal dynamics: 
Resolving chaos 
Two planetesimals in a swarm 


Figure 4 
Figure 5 


Figure 1, Kokubo & Ida (1995) 

Figure B3, Weidenschilling et al. (1997) 


M f ■ 7 ;(/ (: J' a Ig orifJi m : 
Planetary accretion {R = 10^ km) 

{face is the accretion fraction) 
Planetary accretion (i? = 5200 km) 

high-speed collisions 


= 0.138 ±0.002 
face = 0.008 ± 0.001 
Figure 6 


face = 0.140, Grccnzwcig & Lissauer (1990) 
f^^^ = 0.10 ± 0.2, Duncan, Levison & Lee (1998) 
face = 0.009, Greenzweig & Lissauer (1990) 

= 0.008 ± 0.002, Duncan, Levison & Lee (1998) 


Hybrid code: 
Planet formation (coagulation) 
Terrestrial planet formation 


Figures 7-9 
Table 1 


Figures 8 & 9, Weidenschilling et al. (1997) 
Table 1, Chambers (2001) 
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Fig. 1. — The eccentricity of the major planets derive from a symplectic code and the adaptive 
code. With a large timestep dt of one thirtieth of the period of Jupiter's orbit (Pj), the symplectic 
integrator (panel a) generates phase errors — slow drifts from the true phase-space position of the 
planets. These errors are reflected in the relatively limited range of eccentricities, as compared 
with a run with smaller timesteps (panel b). The adaptive code, with a large error tolerance, 
will run with large timesteps. The adaptive results (panel c) generate more realistic eccentricities 
compared to the small-timestep symplectic code and the small-error tolerance adaptive code (panel 
d). However, the adaptive code suffers from a slow drift in energy. The net fractional energy errors 
in the adaptive runs are 10~^ {10~^) for large (small) error tolerances. 
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Fig. 2. — The instantaneous Keplerian semimajor axes of the major planets, with masses scaled 
by a factor of 50, as a function of time from a simulation with our adaptive integrator. Saturn is 
ejected just after 2,000 yr; Neptune has acquired an orbit with eccentricity of 0.92 and a semimajor 
axis of about 250 AU. The overall dynamics are highly chaotic and quite sensitive to the outcome 
of the initial interaction between Jupiter and Saturn in the first 100 yr of the simulation. 
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Fig. 3. — Gravitationally bound binaries orbiting a solar-mass star. The upper panel shows the 
semimajor axis of a pair of Jupiter-mass planets as calculated in the center-of-mass frame of the 
pair as they orbit the star. The binary eccentricity is e = 0.6, as in Duncan, Levison &; Lee (1998). 
The lower panel is the radial position of a test-particle ( "spy satellite" ) on a polar orbit about the 
center of the Earth. In this case, we assume the Earth is spherically symmetric. Each panel shows 
only a single quantity; the appearance of multiple curves is an artifact of sampling. 
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Fig. 4. — Integration of planetesimal orbits. The upper two panels show the orbital separation of 
a pair of 10^^ g objects in units of the mutual Hill radius. In both cases the particles were placed 
on circular orbits near 1 AU with separation of 2.1^/3rH (top) and 1.9\/^rH (middle), where th is 
the mutual Hill radius. When the orbital separation is less than 2\/^rH, then the orbits will cross. 
The lower panel shows interactions between eight particles separated by 2.\^/2>rH- 
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Fig. 5. — Evolution of (e^+i^)^/^ in units of the Hill radius for two lO'^^ g objects and a swarm of 800 
2 X 10^^ g planetesimals. The horizontal axis gives the orbital radius relative to ao = 1 AU around 
a solar-type star. The large filled circles correspond to a pair of the heavy objects at 200 yr (lower 
panel), 800 yr (middle panel) and 2000 yr (top panel). The small filled circles indicate values of 
(e^+i^)^/^ for the planetesimals averaged in radial bins. The x symbols and histograms correspond 
to the two massive objects and planetesimals in a coagulation model (Kenyon & Bromley 2001). 
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Fig. 6. — The minimum mimbcr of timesteps leading to sueeessful pairwise collisions as a function 
of the planetesimal radius. Specifically, nt95 is the minimum number of low-resolution timesteps 
required to resolve a merger between a point-like projectile and a target planetesimal of radius r, 
with a success rate of 95%. The target and projectile are on counterrotating, circular orbits at 
1 AU, initially on opposite sides of the Sun. The error tolerance for the adaptive code is 10~^^ 
(filled circles) and 10~^^ (open triangles). The precipitous drop in n^gs at large r in both cases 
results from the adaptive code's use of high-resolution timesteps as it attempts to account for the 
gravitational effect of the target on the projectile. 
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Fig. 7. — Evolution of particle number (left panels) and eccentricity (right panels) for planetesimals 
orbiting a 1 Mq star at 1 AU. The 32 annuli in the grid contain planetesimals with initial radii 

of 1-4 km, initial eccentricity cq = 1 X 10~^, and a total mass of 0.6 M^. The panels show the 
evolution of N{m) and e for calculations with the coagulation code at t = 10^ yr (upper panels), 
t = 10^ yr (middle panels), and t = 10^ yr (lower panels). Style adapted from Weidenschilling et 
al. (1997). 
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Fig. 8. — Evolution of particle number (left panels) and eccentricity (right panels) for planetesimals 
in orbit around a 1 Mq star at 1 AU, as in Figure 7, but for calculations with the hybrid code. 
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Fig. 9. — Evolution of eccentricity as a function of mass for hybrid calculations at 1 AU. The plots 
show the eccentricity at three different times, as indicated by the legend. For smaller mass values 
(m < 10^^ g), averages are shown for three different annular regions within the computational 
grid. Typically, the eccentricities are higher for objects at smaller annuli, an effect which is more 
pronounced at earlier times. For larger mass values (m > 10^^ g), the eccentricities of single objects 
are shown. These results illustrate the increase in viscous stirring of smaller mass objects with time, 
as well as the dynamical friction generally experienced by the most massive bodies. 
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Fig. 10. — Evolution of the mass-weighted eccentricity e for a hybrid model. The solid circles show 
the evolution of the oligarchs; the open diamonds indicate the evolution of all objects in the grid. 
Viscous stirring and dynamical friction lead to a slow increase in e. Strong dynamical interactions 
produce abrupt rises in e; mergers produce abrupt declines. 
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Fig. 11. — Final eccentricity as a function of final mass for calculations at 0.4-2 AU. The open 
diamonds indicate current values for the terrestrial planets, Mercury, Venus, Earth, and Mars. The 
filled circles plot results for the ensemble of hybrid calculations listed in Table 1. 



